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Abstract: The first processing stage in computational vision, also 
called early vision, consists in decoding 2D images in terms of properties 
of 3D surfaces. Early vision includes problems such as the recovery of mo- 
tion and optical flow, shape from shading, surface interpolation and edge 
detection. These are inverse problems, v^rhich are often ill-posed or ill- 
conditioned. We review here the relevant mathematical results on ill-posed 
and ill-conditioned problems and introduce the formal aspects of regulariza- 
tion theory in the linear and non-hnear case. More general stochastic reg- 
ularization methods are also introduced. Specific topics in early vision and 
their regularization are then analyzed rigorously, characterizing existence, 
uniqueness and stabiHty of solutions. 
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Introduction 

Vision systems, either artificial or biological, are confronted with the 
problem of inferring geometrical and physical properties of surfaces aromid 
the viewer. The available data -the images - consist of two dimensional ma- 
trices of light intensity values measured by an eye or a camera. For tasks such 
as navigation, manipulation and visual recognition, vision systems have to 
recover 3D properties of surfaces from the 2D images. Typical 3D properties 
are the distance between the surfaces and the viewer, their orientation, struc- 
ture, teicture, reflectance and motion parameters (from a temporal sequence 
of images ) . 

The visual skills that provide us with this kind of information have been 
explored in animals and humans with physiological and behavioural tech- 
niques. With the recent development of computer vision, these problems 
have been formulated rigorously and given by now famihar names, such as 
structure from stereo, structure from motion, structure from texture, shape 
from shading, edge detection, visual interpolation, computation of optical 
flow. The computational modules that solve them constitute together the 
core of early vision, and provide spatial and geometrical information about 
the 3D world. The results of this first stage of processing are then used for 
higher level tasks such as navigation in the environment, manipulation of 
objects and of course object recognition and also reasoning about objects. 
Unlike high level vision, early vision is mostly considered as a bottom-up set 
of processes that do not rely upon specific high-level information about the 
scene to be analysed. It is commonly argued on the basis of computational 
and psychophysical considerations that these different modules of early vision 
can be analysed independently of each other, at least to a first approxima- 
tion. Their most natural implementation is in terms of distinct pieces of 
hardware, whose outputs will be integrated at a later stage, possibly using 
more "intelligent" procedures (another possibihty is to use coupled Markov 
Random Field models for the integration stage, see [68]). 

Even a superficial analysis of these problems reveals their common in- 
verse nature: they can be regarded as inverse optics since they attempt to 
recover 3D properties of surfaces from the 2D images they generate. This 
observation characterizes the field of early vision as the solution of problems 
of inverse optics^^'. 

It is well known that inverse problems are very often ill-posedf^i''^i in 
the original sense of Hadamardt"*i ; that is, the solution may not exist or it is 
not unique or does not depend continuously on data. These problems can 
be formulated as discrete or continuous problems according to the type of 



the available data. In order to .so/ye ill-posed problems a priori information 
about generic properties of the solution must be used. Regularization theory 
is a set of techniques that have been developed for this reason. Standard 
regularization exploits a priori knowledge by restricting the functional space 
to which fhe "'solution" belongs: the specific techniques use generalized in- 
verses or variational formulations. An alternative possibihty, that we call 
stochastic regularization, is based on a Bayesian approach and optimal esti- 
mation. All these reasons we introduced recently regularization techniques 
into computer vision. B. Horn (see [69] for a comprehensive review of his 
work) had approached several problems in vision from a similar point of 
view, using minimization techniques for their solution, without an expHcit 
connection with regularization techniques. 

The goal of this paper is to review the mathematical aspects of this 
approach (for a less rigorous discussion, see [Ij). It is organized in two 
parts: Part One reviews the main mathematical results that characterize the 
difference between well-posed and ill-posed problems (Section 2) and between 
well-conditioned and ill-conditioned problems (Section 4). The notions of 
generahzed inverses (Section 3) and of regularization methods (Section 5) 
are then introduced. Section 6 contains some results related to inverse non- 
linear problems. Section 7 illustrates stochastic regularization. 

Part Two shows that several variational principles recently introduced 
in early vision can be formulated as regularized solutions to ill-posed inverse 
problems. Four problems in early vision are studied in detail: edge detection 
and numerical differentiation (Section 8), optical flow (Section 9), surface 
interpolation (Section 10), and shape from shading (Section 11). 



Part One 



1. Outline 



In this part of the paper we review some of the methods which have been 
developed for the approximate sohition of ill-posed problems. The hnear 
case is discussed in detail since a well-developed theory is available. We also 
make some comments on non-linear problems. 

In Section 2 we define the class of well-posed problems, stressing that 
a well-posed problem is not necessarily robust against noise. A well-posed 
problem, in order to have solutions that are robust against noise, must also 
be well-conditioned (see Section 4). For ill-posed, Unear, inverse problems, 
well-posedness can be restored by generaUzed solutions if the range of the 
operator (which has to be inverted) is closed (see Section 3). When the 
range of the operator is not closed, or when the problem is seriously ill- 
conditioned, regularization techniques have to be used (Section 5) in order 
to avoid the instabihty of the solution against noise. Therefore, since images 
are intrinsically noisy, these techniques represent the ideal tool for early 
vision problems. Some results on inverse nonlinear problems are presented 
in Section 6. A stochastic approach to inverse problems is introduced in 
Section 7 and its connections with standard regularization are discussed. 



2. Well-Posed and Ill-posed Problems 

Hadamardt'^i't^] defined a mathematical problem to be well-posed when: 

(a) for each data g in a given class of functions Y there exists a solution 
u in a prescribed class X (existence); 

(b) the solution u is unique in X (uniqueness); 

(c) the dependence of a upon g is continuous, i.e., when the error on 
the data g tends to zero, the induced error on the solution u tends also to 
zero (continuity). 

The requirement of continuity is related to the requirement of stability 
or robustness of the solution (see, for instance, [6]). Continuity, however, is 



a necessary but not siitHcient condition for stability. A well- posed problem 
can be ill-conditioned (see Section 4). 

All the classical problems of mathematical physics, such as the Dirichlet 
problem for elHptic equations, the forward pro!)lem for the heat equation, and 
the Cauchy problem for hyperbohc equations, are well-posed in the sense of 
Hadaniard. Also, the ''direct" problem in scattering (or imaging) theory, 
namely the computation of the scattered radiation (image) from a known 
constitution of the sources and of the targets, is well-posed. 

"Inverse'' problems usually are not well-posed. In most cases an "in- 
verse" problem can be obtained from the '^hrect" one by exchanging the 
role of solution and data. For instance, in the case of scattering theory, 
the inverse problem consists in the computation of the characteristics of the 
targets from the knowledge of the sources and of the scattered radiation. 

Consider a very simple example taken from classical optics. If the energy 
distribution u is given in the object plane of an optical instrument and if the 
characteristics of the instruments are known, it is possible to compute, by 
solving the wave equation, the energy distribution g in the image plane. In 
the case of Fourier optics, one finds a linear relation between u and g: 



g{x) = / K{x,y)u{y)dy, (2.1 



the kernel K{x,y) being the impulse response (point spread function) of the 
instrument. The direct problem (the computation of g given a) is clearly 
well-posed. The inverse problem (the computation of u given g) usually is 
not. 

Assume that K(x,y) = A"(.r -//). where /v(.r) is a band-hmited function. 
Then there exist functions a which produce a zero image (think of a function 
which has only Fourier components out of the band of the instrument) and 
therefore uniqueness does not hold. Furthermore, if g{x) is the result of 
an experiment and therefore is affected by noise, it is not necessarily band- 
hmited or it can have a band broader than that of the instrument. Under 
these circumstances the solution of Equation (2.1) does not exist. 

The need to investigate problems which are not well-posed but are of 
interest in apphed science originated two interesting branches of mathemat- 
ical analysis: the first is the theory of generalized ijiverses f"'-f^i, which 
is an extension of the theory of the Moore- Penrose inverse of a matrix; 
the second is the regularization theory of ill-posed (or improperly posed) 
problems^2l'f^^[*^]'[iol'[i^]. These days, the term ill-posed is used generally 
(but not only) for those problems which do not satisfy the requirement of 



continuity. Examples o{ ill-posed problems are analytic continuation, the 
Cauchy problem for elliptic equations, backsolving the heat equations, su- 
perresolution, computer tomography, Fredholm integral et^uations of the first 
kind, and as we will see, many problems in early vision. 



3. Generalized Inverses 



Most linear inverse problems can be formulated as follows: assume that 
functional spaces X,Y (for instance, Hilbert spaces) are given and that a 
linear, continuous operator L from X into 1' is also given; then the problem 
is to find, for some prescribed g eY, a function u e X such that 

9 = Lu. (3.1) 

In this formulation, the direct problem is just the computation of g, given 
II. Therefore, continuity of L is equivalent to well-posedness of the direct 
problem. Notice that Equation (2.1) is a special case of Equation (3.1). 

The problem (3.1) is well-posed if and only if the operator L is injective 
(i.e.. the equation La = has only the trivial solution u = (uniqueness)), 
and it is onto Y (existence). Then general theorems of functional analysis 
(for instance, the ''closed graph theorem") assure that the inverse mapping 
Z~Ms also continuous (continuity). 

Assume now that the equation Itt = has nontrivial solutions. The set 
of these solutions is a closed subspace of X, which is called the null space 
N{L) of I. Tliis is the subspace of the "invisible objects", since they produce 
a zero image g. Assume also that the range R{L) of L, namely the set of the 
g which are images of some a (E A^ is a closed subspace of 1'. An example is 
provided by the integral operator corresponding to the perfect low pass filter 

IT M \ f ^ smQ{x — y) 

{Lu){x)= — -^ ^u{y)dy. (3.2) 

In such a case, if we take .Y - Y = I'^i -oc, +oo), the null space is 
the set of all the functions u whose Fourier transform is zero on the band 
[-Q,fi], while the range of L is the set of the band-hmited functions with 
bandwidth fi, which is a closed subspace of L^{-oo,+oo). Notice that L is 
a projection operator, the so-called band-Hmiting operator. 

A way of restoring existence and uniqueness of the solution under the 
conditions above is to redefine both the solution space X and the data space 
Y. We take a new space X' which is the set of all the functions orthogonal 
to N{L) (in the case of (3.2), X' is the space of the square integrable Q- 
bandhmited functions), and we take R{L) as the new data space Y' (in the 
case (3.2) again, the space of the square integrable Q-bandHmited functions). 
Then for any g e Y' there exists a unique u t X' such that g -= Lu, (in the 



case of (3.2) the solution is trivial: a = g) and therefore the new problem is 
well-posed. 

The redefinition of the spaces X,Y outHned above usually is quite diffi- 
cult (almost impossible) in practical problems. Therefore, it is useful to have 
a method, based on the solution of variational problems, wliich produces the 
same result. This is just the method of generalized inverses . 

3.1. Least squares solutions or pseudosolutions 

Consider first the case in which L is injective but not onto (i.e., the exis- 
tence condition is not satisfied). The set of functions u e X that solve the 
variational problem 

\\Lu — gWy- = minimum, (3.3) 

where || • \y denotes the norm of F, are called the least squares solutions (or 
pseudosolutions) of Problem (3.1). These solutions can be easily obtained 
considering the first variation of the functional (3.3), 

2Re{Lu - g,Lh)Y, (3.4) 

where h is an arbitrary function of X and (•,-)^ the inner product of the 
Hilbert space Y. Setting (3.4) equal to zero, we obtain the Euler equation 

VLu = L*g, (3.5) 

where L* is the adjoint of the operator I (I* is a mapping from Y into X). 
When R{L) is closed. Equation (3.5) always has solutions but the solution 
is not unique when N(L) is nontrivial. Notice that the set of solutions of 
Equation (3.5) coincides with the set of solutions of the equation 

^« = Pg^ (3.6) 

where P is the projection onto R(L). Therefore, solving Equation (3.5) 
is equivalent to taking 1'' = R{L) or to projecting g onto Y' . When the 
operator L is injective, the solution of (3.5) is unique and well-posedness has 
been restored. 



3.2. Normal pseudosolutions or generalized solutions 

Consider now the case in which L is not injective (i.e., the uniqueness condi- 
tion is not satisfied and the problem is underconstrained). Then, one looks 
for the solution of (3.5) which has minimal norm 

II It II X = minimum. (3.7) 
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This solution is unique and is denoted by u^ . u+ is usually called the gener- 
alized solution (or normal pseudosolution) of Problem (3.1). u^ is orthogonal 
to N{L) and therefore this procedure is equivalent to taking .Y' = A'(L)^. 

Since there exists a unique a+ for any g G F, a linear mapping 1+ from 
1' into -Y is defined by 

"^ = ^^9- (3.8) 

The operator X+ is the generalized inverse of L and it is continu- 
ous. Therefore, the problem of computing the generalized solution of Equa- 
tion (3.1) is well-posed if and only if R(L) is closed. The essential reason for 
this result is that in this case the space Y can be decomposed as 

Y = R(L)®R^(L), (3.9) 

where means direct sum and R^(L) is the orthogonal complement oi R{L). 
This decomposition can be made if and only if R{L) is closed. 

3.3. C-generalized solutions 

In several inverse problems, the generahzed solution is trivial or does not sat- 
isfy some physical requirements such as smoothness. Examples are provided 
in Section 9. Then an extension of the generaHzed solution goes as follows: 
let p(u) be a norm or a seminorm on X (see Appendix A for the definition) 
of the following style: 

p(u) = \\Cu\\^ (3.10) 

where C is a linear operator from X into the Hilbert space Z (the constraint 
space). The operator C may not be defined everywhere on A". For instance, 
suppose X is a space of square-integrable functions and C is a differential 
operator. Therefore, in general, p(u) is defined on a subset of A, i.e. the 
domain of C, denoted as D(C). When the null space of C is trivial (contains 
only the null element of X), then p{u) is a norm on D{C); otherwise, p{u) 
is a seminorm. 

If there exists a unique least-squares solution which minimizes p{u), 
we denote it by aj and we call it a C-generalized solution. The mapping 
g ^ Uc defines a linear operator L^ from Y into A", which will be called 
the C-generalized inverse of I. It is obvious that a^, can have a nonzero 
component onto N(L), the subspace of the "objects" which are "invisible" 
under the action of the operator L. Therefore, this procedure is physically 
plausible only when the constraint describes some physical property of the 
solution of the problem. 
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xN'ecessary and sufficient conditions for the existence of u^ for any g have 
been given in the case where R{L) is closed and C is a bounded operator 
with R(C) also closed^ '^. However, the assumption of a bounded constraint 
operator C may not cover the interesting case of a differential operator. 
Furthermore, when D{C) is a subset of A", it is obvious that a^. does not 
exist for any g t 1'. If we denote by LD[C) the set of all the functions g kl Y 
such that g ^ Lu with u G D{C), then LD{C), in general, does not coincide 
with R{L). Under these circumstances, if Pg ^ LD[C), the intersection 
between the set of the least squares solutions and D{C) is empty and ut 
does not exist. In other words, the problem of determining the C-generalized 
solution may be ill-posed even wdien R{L) is closed. 

Sufficient conditions which assure the existence of » J for any g such 
that Pg e LD{C) are (see Appendix B): 

(i) The intersection oi N(L) and N{C) contains only the null element of 
X, i.e., the set of equations 

La = 0, Cu = (3.11) 

has only the common trivial solution a ~ (uniqueness condition); 
(ii) The operator C : A ^ Z is closed with D{C) dense in X and R{C) = Z; 

(iii) The set of functions u such that g = Lu and Cu = 0, i.e., the set LN{C), 
is closed in Y. 

The third condition is always satisfied in the case of seminorms defined 
in terms of differential operators because in that case N{C) is a finite di- 
mensional subspace of A" and L is a continuous operator. 

When the constraint operator C satisfies conditions (i) - (iii) and fur- 
thermore is bounded, » J exists for any g ^Y and the C-generafized inverse 
L(. is bounded. These properties hold true, for instance, in the case of inter- 
polation problems in reproducing kernel Hilbert spaces (see Appendix C). 

3.4. Generalized solutions for problems with discrete data 

We conclude this section by noticing that problems with discrete data can 
be formulated as (3.1), g being now a n-dimensional vector in an Euchdean 
space. In fact, ignoring the errors in the data, a finear inverse problem with 
discrete data can be formulated as followsL^^I : 

Given a set {Fj"=i of finear functionals defined on A and a set {g,}"^i 
of numbers, find a function u G X such that 

g, ^ F,{u);i =- l,...,/i. (3.12) 
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In particular, when the functionals Fi are continuous on A" by Riesz 
Theorem (see Appendix A), there exist functions il^i, 02,-.-, i\i such that 

F,{u) =(ft,t\)^, (3.13) 

where (•, •)x is the inner product of A'. 

For example, the problem discussed in Section 2 takes this form when 
g{x) is measured in a finite set of points only, say .ri, .r2, ..., .r„, and A" is an 
L space. In such a case we have 

0,(i') = K{i-,,y). (3.14) 

It is important to recognize also that interpolation problems take this form 
when X is a reproducing kernel Hilbert space (see Appendix C). 

This problem is a special case of the problem 3.1 if we consider the data 
gi as the components of a vector g in an N-dimensional Euchdean space 1' 
and if we define an operator L from A" into Y by means of the relation 

[La), = {u,il\)^;i = l,...,r7. (3.15) 

The operator L is not injective: N{L) is the infinite dimensional closed 
subspace of all the functions u orthogonal to the subspace spanned by the 
functions Vi- On the other hand, the range of X, R{L), is closed: R(L) is just 
Y when the functions 0, are hnearly independent, otherwise it is a subspace 
with dimension n' < n. 

Along the hues described above one can introduce generahzed solutions 
or C-generahzed solutions for problems with discrete data. Their deter- 
mination is always a well-posed problem in the strict mathematical sense. 
However, numerical stabihty cannot be guaranteed (see the next section). 

As a final remark, w^e point out that the problem of interpolation by 
means of spline functions can be formulated as a problem of determining 
a generalized or C-generaHzed solution in a suitable Hilbert space (see, for 
instance, [I2j,il3j). A simple example is discussed in Appendix C. 
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4. Well-conditioned and ill-conditioned problems 

As already remarked in previous sections, continuous dependence of the so- 
lution on the data does not yet mean that the solution is robust against 
noise. Generalized solutions of inverse problems with discrete data can pro- 
vide striking evidence of this fact. Therefore, it is necessary to investigate 
more carefully error propagation from the data to the solution when solving 
problem 3.1. 

We assume, as in Section 3, that R(L) is closed, so that the generaHzed 
inverse L^ is continuous. We denote by Ag a variation of the data g and by 
Au^ the corresponding variation on the generalized solution a + . Then the 
standard analysis of error propagation proceeds as follows: 

From Equation 3.8, because of the linearity of £^, we get Au^ = L^Ag, 
which implies 

l|A?i^|[^ ::: |jL"^iliJA5r||^-, (4.1) 

where \\L^\\ denotes the norm of the continuous (bounded) operator L^ (see 
Appendix A). Analogously, from Equation 3.1, with u = {i+, it follows that 

ll^llr < II^H-|I«'"ILy- (4.2) 

Combining Equations 4.1 and 4.2 we obtain the inequaUty 

II Au"*" II V I M II AoIIt^ 

" 11'^ <||i:||||i+||L5Mir. (4.3) 



II^^ILy II^- 

It is important to point out that this inequahty is precise in a certain sense. 
When L is an N x M matrix or L corresponds to an inverse problem with 
discrete data, then equaHty can hold. If L is an operator on infinite dimen- 
sional spaces, then one can always prove that the l.h.s. of Equation 4.3 can 
be arbitrarily close to the r.h.s. 

The quantity 

a = IIIIIIIL^II > 1 (4.4) 

is called the condition number oi the problem. When a is not far from 1, the 
problem is said to be well-conditioned, while when a is large the problem is 
said to be ill-conditioned. 

It is obvious that these definitions are not as precise as that of well- 
posedness. However, what is important in practice is the estimation of the 
condition number since it provides insight into the numerical stabihty of the 
problem. In the case where L is an N x M matrix, ||i:|| is the square root of 
the maximum eigenvalue of the M x M positive semi-definite and symmetric 
matrix L" L (notice that the positive eigenvalues of this matrix coincide with 
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the positive eigenvalues of the matrix LL* ) and ||Z^|| is the inverse of the 
square root of the minimum positive eigenvalue of the same matrix, i.e., 



'mm 

M 



, ^max 



ore generally, if we indicate hy a^{L*L) the positive part of the spectrum 
of the operator X*L, we have 



/majc<7+(L*X) 

In order to provide an example of a well-posed problem which can be 
extremely ill-conditioned, we consider the finite moment problem, i.e., the 
problem of determining a function u[x), defined for example on [0,1], given 
its moment up to the order iV - 1: 



Jo 



Qn ^ I .C" 'u{x)dx\ R = l,...,iV. (4.7) 

Jq 

If we take X = 1^(0,1), then it is easy to recognize^^^l that the operator 
LL* is just the Hilbert matrix 

[HN(-l)\nm = — ■ -; n,m = !,...,.¥, (4.8) 

which is a classical example of an ill-conditioned matrix. From well-known 
results it follows that the condition number for the generaUzed solution of 
problem 4.7 is approximately given by a ^ €xp(l.1^N) and therefore it 
grows exponentially with .V. Already for moderate values of N, a takes 
unacceptable values. 
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5. Regularization Methods 



When the range of I, R{L) is not closed, then the inverse L~^ or the gen- 
erahzed inverse I^ is not defined everywhere on Y and it is not continuous. 
Therefore, both the requirements of existence and continuity do not hold 
true. This is the most difficult case and appropriate techniques are required. 
An example of operators in this class is provided by compact operators (not 
of finite rank) as shown in Appendix A. It is easy to see that an ill-posed 
problem has a condition number a = oo. Therefore, extremely ill-conditioned 
problems behave in practice as ill-posed problems and have to be treated by 
the same technique. 

5.1. Tikhonov regularization 

The most investigated approach to ill-posed problems is the regularization 
method of Tikhonov^^^J. The key idea is to introduce a family of continuous 
"approximations" of a noncontinuous operator. More precisely, a regular- 
ization algorithm for the generahzed solution of Equation (3.1) is given in 
terms of a one-parameter family of continuous operators Rx, A > 0, from Y 
into X, such that for any given g ^ R{L), 

liniRxg = L^g. (5.1) 

Therefore, when apphed to noise-free data g, Rx provides an approximation 
of u^ which becomes better and better as A — > 0. However, when Rx is 
applied to noisy data g^ = g ^ n^ and n^ represents experimental errors or 
noise, we have 

R\9e = R\g + R\n,, (5.2) 

and the second term typically is divergent when A -^ 0. It follows that a 
compromise between "approximation" (the first term) and "error propaga- 
tion" (the second term) is required. This is the problem of the "optimal 
choice" of the regularization parameter A. 

One of the most studied regularization algorithms is obtained by mini- 
mizing the functional 

\\Lu - g\\ y + A II C?/ II 2- = minimum, (5.3) 

where C is a ccnistraint operator, satisfying for instance the conditions stated 
in Section 3. In the original paper of Tikhonov, it is given by 



r=0 
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where the weights Cr{x) are strictly positive functions and il^'^\x) indicates 
the r*''-order derivative of u{x). If a,^ is the solution of (5.3), and if we put 

ux -= Rxg, (5.5) 

then 

i?A = (I*Z-f AC*C)-^I*. (5.6) 

Notice that ux is unique when the Equations (3.11) have only the trivial 
solution a = and that when X -> 0. g e R{L), Rxg converges to Lc^ g '^i. 

Three methods have been proposed for the choice of A in Equation (5.G) 
and in the case of noisy data gc-. 

(i) Among all u such that iCwIiz i: E find u that minimizes \\Lu -geWv ^^^1. 
Using the method of Lagrange multiphers the solution of this problem 
can be reduced to the solution of Equation (5.3), with A arbitrary, and 
to the search of the unique A such that 

il<^^u!lz -= E. (5.7) 

(ii) Among all u such that \\Lu - g,^\Y < £, with given £. find u that min- 
imizes ^\Cu\\z fi^J'Ci^J. Again, the solution of the problem is equivalent 
to finding the unique A such that 

ll^"A ' qAW = s. (5.8) 

This is also called Morozov's discrepancy principle. 

(iii) Among all u such that \\Lii - g^Wv S t, \\Cu\\z ^ E, find a u of the 
type (5.5). This is equivalent f^^^t^^l to taking 

^-(-V^)'. (5.9) 



The first method consists of finding the function u that satisfies the 
constraint \\Cu\\z < E and best approximates the data. The second method 
conputes the function u that is sufficiently close to the data (5 depends on the 
estimate of the errors) and is most ^'regular". In the third method, one looks 
for a compromise between the degree of regularization and the closeness of 
the solution to the data. 



5.2. Regularization and filtering 

The regularized solution (5.5), (5.6) takes a very simple form in the case 
where L is compact and C ^ I the identity operator in X. Then, using the 
singular value decomposition of I (see Appendix A) we obtain 
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2 , . —{9^v^.)yu^., (5.10) 

and therefore the regularized solution is essentially a "filtered" version of the 
non-regularized (generaHzed) solution of Equation (3.1), 

"^ = y\ — (<7,i'fc)F«fc. (5.11) 

k ^ 

This remark suggests that, in this case, one can define regularization algo- 
rithms in terms of filter functions $it.(A): 

ux = y2^^(\) — {g^v^,)yuf^ (5.12) 

— t>t 

k ^ 

satisfying the conditions: (i')*it(A) < 1; (ii')$;t(A) -> 1, for any k, when 

A -> 0; (iii')$fc— bounded for any k and any A > 0. Such a procedure is 

often used in the inversion of compact operators as well as in the inversion 

of convolution operatorsf'^](see Section 8). 



5.3. Smoothing and interpolation 

As already remarked, regularization algorithms can be used also for ill- 
ditioned problems. A well-known example is the smoothing of a function 
whose values, specified on a finite set of points, are affected by errorsf^^'. It is 
interesting to compare smoothing and interpolation by means of cubic splines 
using the framework outlined above. Interpolation of a function {/(.r),.r 6 
[0,1], is the problem of searching for a function which takes the prescribed 
values 

ii{'P,)=gi ; i = l,...,n (5.13) 

and minimizes the seminorni'^"^' 

p{u) - / \u"{x)\'^dx. (5.14) 

Jo 

Therefore, the interpolation problem is equivalent to the computation of a 
generalized solution. On the other hand, the smoothing problem is formu- 
lated again as the minimization of the seminorm (5.14), but condition (5.13) 
is replaced by 



Y^H'^i. 



9i\^ < £^ (5.15; 



i=\ 
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(for simplicity, we have assumed that the errors on the data have the same 
variance). Therefore, the smoothing problem corresponds to method (ii) for 
the choice of the regidarization parameter. 

5.4. Cross validation and generalized cross validation 

We conclude this section with a short description of the cross validation 
mef/io</f*^^l'f^'^l . This is a method for the choice of the regularization parameter 
and it has been applied to smoothing problems and also to the solution of 
Fredholm integral equations of the first kind in the framework of the method 
of collocation (or moment-discretization). However, it appHes to any hnear 
inverse problem with discrete data, as formulated in Section 3. 

The idea behind cross vahdation is to allow the data points themselves 
to choose the value of the regularization parameter by requiring that a good 
value of the parameter should predict missing data points. In this way, no a 
priori knowledge about the solution and/or the noise is required. 

Let {Lu)i be defined as in (3.15) and let u^^^ be the minimizer of the 
functional 

F^xh^A = lf^=l\{Lu),-g,\'+\\\u\\],. (5.16) 

Then the cross vahdation function Fo(A) is defined by 

^'^(^) = ^EK^"a'U-^^-I' (5.17) 

fc=i 

and the cross validation method consists in determining the value of A, say 

A, which minimizes (5.17). The computation of the minimum is based on 

the relation 

where ux is the minimizer of the functional 



1 " 
F\u\ = -V|(Iw), -i// +A| 

k=\ 



12 



;5.i9 



and Afcfc(A) is the kh-iXx entry of the n x n matrix 

.4(A) ^LL\LV ^\IY\ (5.20) 

where LI* is the Gram matrix of the functions ^ (see Equation 3.15). 

It has been shown '^^^ that, from the point of view of minimizing predic- 
tive mean square error, the minimization of Vo(A) must be replaced by the 
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inininiization of the generalized cross-validation function defined by 

V{\) ^{-Tr[I ^- A(\)r^-\\(I - A{\))g\\'), (5.21) 

ri H 

where i| -jl denotes the EucUdean norm and Tr is the trace operation. An 
important property of V'(A) is the invariance with respect to permutations 
of the data. 
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6. Regularization of nonlinear problems 

The case of nonlinear ill-posed problems is quite difficult and, for the mo- 
ment, no general approach seems to exist. 

If A is a nonlinear operator from a Hilbert space X into a Hilbert space 
1', we have the equation 

9 = A{u). (6.1) 

Obviously, a solution of this equation exists if and only if g is in the range 
of the operator A. 

6.1. Linearization 

The simplest way of treating Equation (6.1) is to try to Hnearize the problem. 
This is the case of a differentiable operator^'^'^^ . The nonhnear operator .4 has 
a first derivative at the point Uo if there exists a Hnear operator Lo : A" — > 1' 
such that, for any u £ X, 

Um-[A{uo + tu) - A{u^)\ = LoU. (6.2) 

The operator Lo is called the first derivative of A at the point Uo and one 
usually writes 

Lo = A'{uo). (6.3) 

An operator which is differentiable at the point Uo is also continuous at that 
point. 

If an approximation Uo of the solution of Equation (6.1) is known and 
if the operator A is differentiable at Uo, then Equation (6.1) can be approx- 
imated by the linear equation 

dgo = Loduo, (6.4) 

where dgo = g - A(uo), diio = u - Uo, and Lo is the derivative of A at Uo- 
Obviously, the procedure is consistent if the solution duo of Equation (6.4) 
is a "small" correction to the approximate solution iio- 

The procedure can be iterated. By means of the solution duo of Equa- 
tion (6.4), one gets a new approximation, iii = Uo + duo, of the true solution 
{/. Then one considers the linear equation dgi = Lidui, where Li = A'{u.i ), 
dgi = g - -4(fii), and dui = a - ui. By solving this equation one gets a new 
approximation ^2 = wi f ^"i and so on. It is easily recognized, by writing 
Equation (6.1) in the form P{u) = with P{u) = A{u)-g, that this method 
is just an extension to functional equations of a method which, in the case 
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of real equations, is known as Newton's method or the method of tangents. 
Such an extension is also known as the Newton-Kantorovich method and it 
is one of the few practical methods for the actual solution of a nonHnear 
functional equation. 

The iterative algorithm can be put in the following form: 

"n + l = Un ^ [A'(Un)]'^[9 - A{Un)]; (6.5) 

and a simplified algorithm is given by 

Un+l = Un + [A'{Uo)]'^lg - A{Un)]. (6.6) 

Sufficient conditions for the convergence of both iterative algorithms have 
been given^'^^ They include the continuity of the inverse of the derivative 
of the operator A. In several inverse problems this condition is not satisfied. 
It has been suggested^"! to use, at each step of the algorithm, a regularized 
approximation of the inverse of the derivative of the operator .4. Convergence 
results for such a modified algorithm are not yet available. 

6.2. Generalized and regularized solutions 

Extensions of regularization theory to ill-posed nonHnear problems have also 
been proposed: the case of nonHnear integral equations has been investigated 
by Tikhonovt^^l and an abstract approach is given by Morozov^^'] . 

We assume that A : X -^ Y is a continuously diff'erentiable operator, 
i.e., that A has a derivative at each point u G -Y and that this derivative is a 
Hnear continuous operator. However, even in the case of such a simpHfying 
assumption, a well-developed theory of generalized inverses does not exist. 
One can introduce least-squares solutions of Equation (6.1) by solving the 
variational problem 

\\A(il) - qIy = minimum, (6.7) 

analogous to the problem (3.3). When a solution of such a problem exists 
for any g G Y, one says that Equation (6.1) is strictly normally solvable. 
A sufficient condition for strict normal solvabiHty is that the range of A is 
weakly closed in Y f^^'. Notice that this condition may be stronger than the 
condition of closure of the range which appHes to the case of Hnear operators 
(Section 3). Weakly closed sets are (strongly) closed, but the converse is not 
always true. 

If, for a given g, the set of least squares solutions is not empty, one could 
try to select one of these solutions by means of another variational principle 
as in Section 3.1; i.e., by minimizing a norm or seminorm such as (3.10). 
In contrast to the case where the operator A is Hnear, the generaHzed or 
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C-generalized solution defined in such a way may not exist and, even if it 
does exist, is not necessarily unique. Such a lack of uniqueness appUes also 
to the case of regularized solutions (in which case, however, existence can be 
easily assured). 

The basic point in the definition of regidarized solutions is again the 
minimization of a functional similar to (5.3); i.e., 

*aW = \\A{u)-g\\l^ + \\\Ca\\l. (6.8) 

The uniqueness of the minimum of $a[«1 usually is not proven (but see [26] 
for a special case where uniqueness holds true). However, it is not difficult 
to prove the existence of at least one local minimum. Here we give the proof 
under conditions which are satisfied in the case of the problem of shape from 
shading (Section 11). 

Assmne that the operator A : X -^ Y is continuous everywhere and 
that the constraint operator C : X -^ Z is linear and has a compact inverse 
C~ . (This condition is satisfied, for instance, by the difi"erential operator 
(5.4)). Then, for any A > 0, the functional (6.8) has at least one minimuin 
point u\. The proof goes as follows: 

Let Ma be the lower bound of ^\[u\ and let {w„} be a minimizing 
sequence such that 

Ma < ^x[un] < Mx + l/n. (6.9) 

It follows that: 

lCu„\\z<{\i^[u„]Y''<(^^Y", (6.10) 

and therefore, the sequence {vn = Curr} is bounded. Since C^Ms compact, 
we can extract from {un = C"^r„} a subsequence strongly convergent and 
such that the corresponding subsequence of the Vn is weakly convergent. 
Without loss of generality, we can assume that these conditions are satisfied 
by the sequence {a„} itself. Then, let ux be the strong limit of {u„} and ^a 
be the weak limit of r„; it follows that ux = C'^vx- 

Since Cun weakly converges ioCux, from the lower weak semicontinuity 
of the norm we have 

WCiLxWz < liminf ||Cu„||^. (6.11) 

On the other hand, from the strong convergence of u^ to ux and from the 
continuity of A{u) we have 

\\A{ux) - gWr = lim ||A(m„) - 5r||v', (6.12) 

and, by combining Equations (6.9), (6.11), and (6.12), we get 
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Ma < ^[u\] ^ limiiif #A[«'nl = lim$A["nl = ^a. (6.13) 

It follows that 

^[ux] = Ma, (6.14) 

and the existence of the minimum point is proven. 

As stated above, in general nothing can be said about the uniqueness of 
the minimum of the functional (6.8). However, if we assume that: 

(a) for a given g, Equation (6.1) has a unique solution {/ in the domain of 

(b) in a neighborhood of u, the operator .4 has everywhere continuous first 
and second derivative; 

(c) the derivative of A at a, A'{u), is invertible; 

then, by a rather easy generaUzation of the theorems contained in [26], one 
can prove that if g^ is noisy data, with '\g ~ geliY- -- "' ^^^^ i*" i^i ^^^^ f^mc- 
tional (6.8), with g replaced by ^f., we choose the regularization parameter 
A in such a way that A = 7t^ where 7 is an arbitrary constant, then any 
minimum point of such a functional converges to a when t i-^ 0; therefore, 
for sufficiently small values of t, there exists only one minimum point. 
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7. Stochastic route to regularization 



When a priori knowledge of statistical properties of the signal and of the noise 
is available, a probabiHstic version of regularization methods is possiblet^^l't^^l' 

[31]_ 

Here we consider a Bayesian approach that has the advantage of showing 
the connection between Markov Random Field models and standard regu- 
larization as developed in this paper. In particular, we will be able to see in 
which sense standard regularization is a special case of MRF models and is 
itself equivalent to Wiener filtering. 

The first step is to write Equation (3.1) in this form, 

g = Lu ^ w, (7.1) 

where w is a function representing the effect of the noise on the data. Notice 
that in this representation no assumption is impHcit about the structure of 
the noise (additive noise, signal dependent noise, etc.). 

The second step is to assume that there exist stochastic processes u,g,w, 
related by 

l=Lu + w, (7.2) 

and that the functions g, u, it? appearing in Equation (7.1) are values of the 
processes u, g, w, associated with a specific outcome of a given experiment 
(we use here the nomenclature introduced in [32]). 

For simphcity, we also assume that the processes u, </, w have zero mean. 
This assumption is in fact not restrictive because, if it is not true, it is always 
possible to introduce processes satisfying this condition just by subtracting 
the means from the original ones. Thanks to the Hnearity of Z, relation (7.2) 
holds true also for the new processes. 

When the mean is zero, the autocorrelation and the autocovariance of 
a process u coincide. If the values of the zero mean process u are functions 
of a variable .r (eventually multi-dimensional), the autocovariance function 
of u is 

Cu_{r,x')^E{u(.r)u{x')}, (7.3) 

where E indicates the expected value. As in the previous sections, we assume 
that the functions u, the values of the process u, belong to a Hilbert space A" 
(for instance, a space of square integrable functions) and that the functions 
g,w, values of the processes g^, iv respectively, belong to the same (possibly 
different from X) Hilbert space Y. The appropriate description of stochastic 
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processes with values in Hilbert spaces is given in terms of weak random 
variables or cylinder set measures f'^'^L 

Then, the autocovariance function of the process u can be considered as 
the kernel of an operator R^ defined on the space X: 

{Ru<?){^v.) = Jc^{x,x')cl>{x')dx\cf> e X. (7.4) 

The operator Ri^ is called the covariance operator (or the covariance) for 
the process u. It can also be defined for weak random variables with values 
in an abstract Hilbert space A" t^^' . 

Coming back now to our basic equations (7.1), (7.2), the inverse problem 
consists in estimating a value of u, given an observed value g of ^ and given 
a priori probabilistic knowledge on the processes u and w. 

We take a Bayesian approach and write the a posteriori probabiHty 
density as 

P{u/g) = constP[a)P{g/u) (7.5) 

where P{u)\s the a priori probabiHty density of process u and P{g/u) is the 
conditional probability density of the data g given u. 

We consider now the special case of u being a gaussian process (or 
equivalently the linear transformation - such as a derivative - of a gaussian 
process). In this case, the a priori probabiHty distribution of jx is 

P{u) = const -expl-- {u,R~ hi) ]. (7.6) 

Let us assume that the noise process w is additive, white and gaussian 
with variance a^ . Then the a priori probabiHty P{g / a) can be written as 

P{glu) =con.s^-exp[--— 11^ - ^^'irv'l- (7.7) 

Depending on the optimality criterion there are now several ways of 
obtaining the best estimate of it given the data g. A commonly used estimate 
is the Maximum A Posteriori (MAP) estimate 

P^'^hesild) = max{P(i//(7)|i/, G A}. (7.8) 

From Equations (7.5) - (7.7) we have 

P(u.lg) ^ 

const ■ ex p[- — (\\Lii - g\\\- + (t\u,R^\l)x)\. (7.9) 

If we put 

R.^^iC'^Crh (7.10) 



25 



then from Equations (7.8) and (7.9) we have 

.U(»|3Pg^) =: mm{A/(u)|« e X}, (7.11) 



where 



M{ii) = \\Lu ~ g\\y \- G-%Ciify^. (7.12) 

It follows that 1%^^^ = F^g where Fq is given by (5.6) with A = cr^ If 
we put R^ = a^I, where I is the identity operator in Y (/?„, is the covariance 
operator of white noise), then Fq can also be written in the following form: 

Fo - RaL*{LR^* + R^~\ (7.13) 

with Ru given by Equation (7.10). 

The operator F^ is sometimes called Wiener filter (or Wiener- Kolmo- 
goroff filter) and is quite similar to the operator (5.6)^^^^. In other words, the 
regularizing operator (5.6) is equivalent to a Wiener filter in the case of white 
noise, provided that the constraint operator C is related to the covariance 
operator /?„ by the relation (7.10). 

Most of our previous assumptions can be relaxed in a more general prob- 
abilistic scheme based on the formalism of Markov Random Fields (MRF) 
defined on finite lattices. In particular, the noise may not be additive, the 
operator L may not be hnear and P{u) does not need to be gaussian. 

MRF models have been formulated for several problems in early vi- 
sion. Under simphfyiug assumptions they reduce to the discrete equiva- 
lent of standard regularization. Though they are computationally expensive 
they may represent a powerful extension of the methods described in this 
paperf^^'l't^^l'f^^]. 
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Part Two 

The initial stage of machine vision, now called early vision, consists of 
distinct but interrelated problems like "edge detection," "computation of 
optical flow," "structure from motion," "structme fnun stereo matching," 
etc. From a theoretical point of view, these problems can be considered as 
independent, at least to a first approximation. The integration of various 
outputs is performed at a higher stage, where geometrical reasoning and 
much a priori information will be used. These different modules may reflect 
processing occurring in our brain, where sinndtaneously we compute different 
information from images: we can extract rapid changes in image brightness 
(edge detection); we can recover the shape of an object from its shading 
(shape from shading); we can understand the motion of objects from the 
changing images (computation of optical flow); we recover the 3D structure 
of a scene from a pair of images (structure from stereo); and we are able to 
have a dense description of 3D surfaces from sparse features (visual surface 
interpolation). 

Several of these problems have been recently solved with variational 
teclmiques, in particidar by Horn, Crimson and Hildreth. We will show that 
many of these results and several new ones - in particular existence and 
uniqueness of solutions - are direct consequences of mathematical results 
presented in Part One. 

Part Two is divided in four sections, each dealing with a main topic of 
early vision. Section 8 presents the ill-posed nature of numerical differenti- 
ation. In Section 9 we show how recently obtained mathematical results on 
optical flow [^'^]'[38lJ39],ji5j^j,p straightforward consequences of regularization 
theory. Section 10 discusses a recent approach to surface interpolation, il- 
lustrating how variational principlesL^'^J'f'^'^J't-^^J't'*^^''-*^]'^'*-*! can be viewed as 
regularized solutions to discrete ill-posed problems. Section 11 reviews recent 
variational approaches to shape from shading, in the framework of regular- 
ization theory. 
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8. Edge detection and numerical differentiation 



Recently standard regularization techniques have been appHed to a classical 
problem of early vision - edge detection. Edge detection, intended as the 
process that attempts to detect and locaHze changes of intensity in the image 
(this definition does not encompass all the meanings of edge detection) is a 
problem of numerical differentiations^^'. As we will see in the next section, 
differentiation is a common operation in early vision and is not restricted to 
edge detection. Differentiation is ill-posed because the solution does not de- 
pend continuously on the data. The intuitive reason for the ill-posed nature 
can be seen by considering a function f{x) perturbed by a very small (in 
L2 norm) "noise" term esinHo^. f{j;) and f(x) + e sin ri.c can be arbitrarily 
close for very small e, but their derivatives may be very different if Q. is large 
enough. This simply means that differentiation "amplifies" high-frequency 
noise. Differentiation can also be seen as the recovery of the solution u to 
the inverse problem g — Lu where 



L:X ^ Y {Lu)(x) = r u(y)dy. 



(8.1) 



Thus u. is the derivative of the data g. The operator (8.1) is not bounded in 
L^( -00, 00), and the range of L is not closed. Therefore, the inverse problem 
is ill-posed . 



8.1. Regularization of different iat ion 

As shown in Section 3, it is possible to restore well-posedness by redefining 
the solution space X and the data space Y. Let us redefine the solution 
space X as the subset X' of square integrable functions f(x) in (-00, +00) 
such that 



-00 



+ 00 , 

[l^—)\F{^)\'duj (8.2) 



exists, where F(u;) is the Fourier transform of /(.c). The new data space 
F' is simply the range of L. It is easy now to see that the inverse problem 
L : A"' — > Y' when the operator L is defined as in Equation 8.1 is well-posed. 

Differentiation can be transformed into a well-posed problem by using 
Tikhonov's regularizing operators. For equations of the convolution type, 
such as (2.1) with K(x,y) = K[x-y), the regularizing operators correspond 
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to convolving g{x) with a filter /(.r, A) (where A -> is the regiilarization pa- 
rameter) whose Fourier transform F(u;,X) satisfies the following conditions: 

(i) < F(u;, A) < 1 for A > and all uj; 

(ii) F(u;, A) is an even function with respect to u; and 6 Lzi -oo, foe); 

(iii) F(u;,A)ja; belongs to l2(-oo, -f oc) for any A > 0; 

(iv) limA^o^(u;,A) = 1. 

This regularizing filter is equivalent to a smooth low pass filter. Three 
main types of filtering have been used in edge detection. We will Hst their 
main properties below. 

8.2. Band- limited filters 

Band-hmited filters are an obvious choice for regularizing differentiation, 
since the simplest way of avoiding harmful noise is to filter out high frequen- 
cies that are amplified by differentiation. Linear and circular prolate func- 
tions constitute an interesting class of band-hmited filtersf'*^''^^^! and have 
already been used in edge detectionf^^^' . These filters satisfy all conditions 
of Tikhonov needed to regularize differentiation, if we take the inverse of the 
band-width as the regularization parameter. 

8.3. Support-limited filters 

All real filters have a finite extension and are support-limited. A class of 
support-Uiiiited filters that has been used in edge detectionf^*^' is the so- 
called difference of boxes (DOB). These filters are Haar functions f^^^ which 
form a basis for square integrable functions on a bounded interval. However, 
these filters do not satisfy Condition (iii) above, and therefore cannot be 
used to regularize differentiation. This conclusion derives from the fact that 
the Haar functions are discontinuous. As a consequence, the hmit of their 
Fourier transform for a; going to infinity tends to zero as uj~^. 

It is possible, however, to introduce smooth support-Hmited filters whose 
Fourier transform tends to as desired as u7 ~> oo. If the function /(.c, A) has, 
for instance, continuous derivative up to order p and the (p+1 )*^^ derivative is 
integrable, then F{uj,X) tends to zero as |u;p^^^^^ Furthermore, if / (.r, A) 
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is C^, then F {u;, A) tends to zero more rapidly than any inverse power of 
u:. An example is provided by the function 

lO, |.r| > A 

where C\ is a constant such that F(0, A) = 1. Therefore, the best support- 
limited filter for edge detection and numerical differentiation is not the DOB 
but the filter (8.3), which is often used in digital signal processing when 
aliasing needs to be reduced. 



8.4. Filters with minimal uncertainty 

The Gaussian function minimizes the product of spread in the space and 
in the frequency domain^^^' and can be viewed as a filter with minimal un- 
certainty. Filtering with a Gaussian function regularizes differentiation, be- 
cause the Gaussian function /(.r, A) = exp(-.cVA^) satisfies all conditions of 
Tikhonov. Moreover, filtering with a Gaussian transforms a continuous and 
bounded function into an entire function. 

8.5. Interpolation and approximation 

Numerical differentiation can also be regularized in a different way. It is 
possible to interpolate or approximate the data with an analytic function 
and subsequently compute the analytical derivative of the interpolating or 
approximating function. 

For instance in ID the "image" model is y, = f{x,)i- e,, where y, is the 
data and e, represent errors in the measurements. We want to estimate /' 
from an interpolating or approximating function /. We can choose a regu- 
larizing functional ||C/f = J{f"{x))^dx, where /" is the second derivative 
of /. This choice corresponds to a constraint of smoothness on the intensity 
profile. Its physical justification is that the (noiseless) image is indeed very 
smooth because of the imaging process: the image is a bandHmited func- 
tion and therefore has bounded derivatives. We can decide that the data 
are noiseless and therefore we want to interpolate the data. This procedure 
is equivalent to interpolating the data with cubic sphnes and differentiation 
can be safely obtained by the analytical differentiation of the interpolating 
spline. If the data are noisy and we want to take errors in the measurements 
into account, we can look for an approximating function minimizing 



E(^^ - /(•^^))' + A / {f'Mfdx. (8.4) 
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In [52] it was shown (a) that the solution /(.r) can he obtained by 
convolving the data y, (assumed on a regular grid and satisfying appropriate 
boundary conditions) with a convolution filter R, and (b) that the filter R 
is a cubic sphne with a shape very close to a Gaussian and a size controlled 
by the regularization parameter A. Differentiation can then be accomphshed 
by convolution of the data with the appropriate derivative of this filter. The 
optimal value of A can be determined for instance by cross vaUdationf^^''^^-^' 
and other techniques. This corresponds to finding the optimal scale of the 



filt 



er 



These results can be directly extended to two dimensions. The resulting 
filters, which are spline filters for discrete data and Butterworth-Uke filters 
for continuous data (they are eventually indistinguishable in practice) are 
very similar to the derivatives-of-a-gaussian extensively used in recent years 

[46], [53], [541, [55] T-f iU „ 1 • • r i.- i lly-»i-ii • 

' ^ ^ ' ^ ^"^ '.It the regularizing functional ||C/|| is 

/ / (V-grad/)%/j- dy, (8.5) 

where V^ indicates the Laplacian and grad/(.r,/y) the gradient of /(.r,^), it 
has been sliownt^"' that the solution f{.i\y) can be obtained by convolving 
the data g(x,y) with the filter 

R2{^.y) = :: / ^-^-—ajduj. (8.6) 

2 Jq Au;*' 4- 1 

where Jq is the zero order Bessel function and z = ^/.r^ r y^ . If the 
regularizating functional is 



the filter becomes 



// 



(V'f(x.i/)fdxdy, 



2 Jo Au;* -f 1 



Therefore numerical differentiation can be regularized in a number of 
ways which are all consequences of the results presented in Part One. There 
are two main possibihties: filtering the data with appropriate derivatives of 
Tikhonov filters; or interpolating (or approximating) the discrete data with 
sphnes and then performing an analytical derivation. These two regidarizing 
procedures are equivalent. 
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9. Computation of optical flow 

A major aim of early vision is the segmentation of the visible scene in regions 
corresponding to distinct rigid objects. Motion is an important source of 
information for this goal. The imaging device projects the 3D velocity field 
of viewed objects into the image as a 2D field. When a moving scene is 
viewed by a camera it is possible to recover directly the optical flow. The 
optical flowi^*^! is conmionly defined as the distribution of apparent velocities 
of movement of brightness patterns in an image. Optical flow and the 2D 
motion field are related and their relationship has been carefully analysed in 
[451. 

In this section we discuss two different approaches to the computation 
of optical flow. Horn and Schunck^^^^l derived equations relating the change 
in image brightness E{.r,y,f) at a point {x,y} and time f to the motion 
of brightness pattern. Their key assumption is that the brightness of a 
particular point in the moving pattern is constant, so that the total derivative 
of E{x, y,t) is zero: 

— (.r,i/,/) = 0. (9.1) 

Then, from local measurements of the partial derivatives oiE(x,y, t) with re- 
spect to space coordinates and time, it is possible to estimate the component 
of the velocity field parallel to the gradient of E{x,y,t). The normal com- 
ponent is not determined and it must be recovered (see [45] for an analysis 
of the validity of the underlying assumptions). 

Hildreth^-^^''^-^*^^ suggested computing the optical flow not over the entire 
image but only along 1-D contours. In real images, these 1-D contours are 
edges corresponchng to sharp changes in image brightness (see Section 8). 
Hildretht^^l''^*' observed that it was possible to obtain the normal vectors 
along the contour by a simple inspection of the extracted edges: if E(x,y, t) 
is again the image brightness, then the normal component v-^ of the local 
velocity vector P" at the points of the contour F is given by 

where V^ is the Laplacian. A better estimate of r^, however, is 



1^ _ d d^E 



where d^ / dn^ is the second derivative along the direction of the gradient^^^^J. 

In this section we will discuss the ill-posed nature of the recovery of 
optical flow as proposed by these authors. 
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9.1. Optical flow along a contour 

We first consider the problem of determining the two-dimensional optical flow 
along a contour T in the image assmning that local motion measurements 
along the contour provide only the component of the velocity in the direction 
perpendicular to the contour. The component of velocity tangential to the 
contour is invisible to local detectors that examine a restricted region of the 
contour. The local velocity vector V{s) is decomposed into a perpendicular 
and a tangential component to the curve 

V(s) = v^{s)t^v^{s)n. (9.4) 

Here s is the arclength and t,n are unit vectors respectively tangent and 
normal to the contour F 

-^ / cos 6* \ . / - sin \ 

'^IsinJ ''=[ cos^ j' (9-^) 

where is the angle between t and the unit vector of the ^--axis. They depend 
also on s but we omit this dependence for simpHcity of notation. 

The component v^is) and the vectors f,n are given by direct measure- 
ments and therefore are the data of the problem. We will denote by g{s) the 
measured values of v^{s) and by g(s) the corresponding velocity field 

9(s)=g{s)n. (9.6) 

Then the problem can be formidated as the inversion of a projection operator 
in the space A^ = Y = L^F) ^^(r) (L'^T) denotes the space of square 
integrable functions defined over T). The norm of a velocity field F G X is 
defined by 

WVWx = J V{s)-V{s)ds = 

J \v^{s)\'ds+ f \v'-(s)\'ds. (9.7) 

The projection operator is 

LV{s) =v^{s)ri, (9.8) 

and the set of the solutions of the equation 

LV = g, (9.9) 

with g defined by Equation (9.6), is the set of the velocity fields V given by 

V(s) ^ilis)t + g{s)n, (9.10) 

where g(s) is the given data function and ir{s) is an arbitrary function in 
X^(r). The generaUzed solution, or solution of minimal norm, exists for any 
data function g{s), but it is trivial since it is given by 

V^ = g (9.11) 
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In other words, the generalized solution restores well-posedness, but it gives 
a solution which does not have any physical relevance. Therefore, one has to 
look for suitable C-generahzed solutions corresponding to physically accept- 
able velocity fields. 

9.2. A C-generalized solution for the optical flow along a 
contour 

A seniinorni introduced by Hildreth gives a very useful constraint for the 
recovery of the optical flow. Put Z - X ^ I''^(r) 9 L-(r) and introduce the 
operator 

CV = V (9.12) 

where the dot means derivation with respect to .s. Then the C-generalized 
solution is the velocity held of the form (9.10) which mininuzes the functional 

IICV'll'^- = / V'-fcLs. (9.13) 

It is easy to show- that existence and uniqueness of the C-generalized 
solution can be derived from the general result given in Section 3.3. 

First, consider the question of uniqueness. We know that the C-gener- 
aHzed solution is unique if and only if the intersection of iV(C) and N[L) is 
the null element (Condition (i) of Section 3.3). Now N[C) is the set of the 
constant velocity fields (or translations), say V = a. Furthermore, iY(L) is 
the set of the velocity fields orthogonal everywhere to /T, i.e., v • n ~ 0. This 
condition can be satisfied by a ^ only if n is constant; that is, only if P is 
a straight fine. Therefore if P is not a straight line, the intersection of N[C) 
and N[L) is always the null element, and uniqueness is restored by the use 
of the C-generalized solution (9.13). 

The existence of the solution follows from the fact that the operator 
(9.12) satisfies Conditions (ii) and (iii) of Section 3.3. Condition (ii) is a 
rather general property of differential operators (see Appendix B for com- 
ments), and Condition (iii) is also verified })ecause N[C) is a two-dimensional 
subspace of X = L'^(P) tp I^(r). Therefore, we can conclude that the C- 
generahzed solution exists whenever </' e LD{C). 

In order to see more precisely the meaning oi the last condition, assume 
that the contour P consists of a finite number of regular arcs, so that the 
tangent is continuous on P with the exception of a finite number of points, 
si , 5-2, ... , .s„, where the tangent has both right and left finiit. Then a solution 
V(s) of the form (9.10) is in the domain of the constraint operator C if 



34 



il'{s)t and g{s)u are differeiitiable on each regular arc and furthermore they 
satisfy suitable conditions at the discontinuity points si in order to assure 
the continuity of V(s). We can derive these conditions from the equations 

V\{s,) = F_(5,); I = l,...,n, (9.14) 

where -\- and - denote respectively right and left Hniit. It follows that 
i'+i^i) = (sin</),)-^ [g+{s,) cos (I), - g-(s,)] 
0-(5,) = {sm(l),)-'^[g+(s,) -g^{s,)cos(p,], (9.15) 

where sin<^ = f_ • 77+ = -f+ ■ 77_,cos 4) = n+ - n_ = f^ • t_ (f^ is the right 
limit of the tangent, etc.). Therefore, if g{s) admits a right and left hmit at 
the points si, it is possible to derive from Equation (9.15) the right and left 
Umit of 0. AH these conditions characterize the subset D{C) which contains 
the unique solution which minimizes the seminorm (9.13). Of course, if ^ is 
not differentiable on the regular arcs or does not have left and right limits at 
the discontinuity points, the C-generalized solution does not exist. It follows 
that the problem is ill-posed. 

Before discussing this point, we want to point out that, if the data 
g are not affected by noise, the C-generalized solution coincides with the 
true solution in two important casest^^^C^^]; the first is a translation of an 
arbitrary contour and the second is an arbitrary motion of a rigid polygon. 
These results can be derived from the Euler equation for the C-generaHzed 
solution. 

Assume that the regular arcs have a differentiable curvature. From the 
following relations, which are true on each regular arc 

f= On , n = -St (9.16) 

where 9 is just the curvature, one can derive from Equation (9.10) 

V(s)= [iis)-e{s)g{s)\t+ [g{s) ^ e{s)iis)]n (9.17) 

and therefore, when V' satisfies the conditions (9.15) 

\\cv\^x = j {\9{s)\' + \e{s)g{s)\'}ds+ 

+- I {\iis)\' + \e{s)ilis)\' + 2e(s)g(s)il;{s) - 2e{s)g{s)iis)}ds (9.18) 

This is a functional of 0, which is an arbitrary function except for Ijeing 
differentiable and satisfying conditions (9.15). Then, by annihilating the first 
variation of this functional, it follows that, on each regular arc, the function 
which minimizes the functional is solution of the differential equation 

-iis) + \e{s)\'iiis) + 2e(s)g(s) + e{s)g(s) = o. (9.19) 
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111 the case of a closed contour, the C-generahzed solution is given by 
the unique solution of Equation (9.19) satisfying the conditions (9.15). If 
the contour is regular everywhere, then one has to add boundary conditions 
such as 

V'(O)-0(/) , V'(0) = ^'(/). (9.20) 

W hen the contour is open, one needs boundary conditions at the end points 
of the contour. These can be obtained directly, through a partial integration, 
from the annihilation of the first variation of (9.18) 

0(0) = ^(0)^(0) , m^O{l)g(l). (9.21) 

However, these conditions are correct only in the case of a pure translation. 
In the general case it is necessary to measure the tangential velocity of the 
end points and take 

0(O) = i'^(O) , 4il) = v^{l), (9.22) 

where r^(0) and v^{l) are the measured values. 

If the motion of the contour is a pure translation and a = {cii, 02} is the 
constant velocity field, the noise-free data are given by 

g{s) = -ill sin ^ 02 cos 0. (9.24) 

Then, if we put 

V'(^) = ai cos^ + 0-2 sin^, (9.25) 

taking into account that ^ = 9g^ g = -O^,^ it is easy to verify that 
satisfies Equation (9.19). In the case of an open contour, also the boundary 
conditions (9.21) are satisfied (the boundary conditions (9.15) are obvious 
since the velocity field is continuous). 

In the case of a rigid polygon, since an arbitrary rigid motion consists 
of a translation plus a rotation, on each segment of the polygon both the 
normal and tangent velocity are linear functions of the arclength s. But, on 
a segment of a straight Hne, Equation (9.19) becomes 4is) = and therefore 
0(5) is a Hnear function of 5. The boundary conditions (9.20) (plus the 
boundary condition (9.22) in the case of an open polygon) give the correct 
values of the constants provided that also in this case the measured values 
are noise-free. 

As we already remarked, the difficulty of this approach is that the prob- 
lem of determining such a C-generaHzed solution is ill-posed. For this rea- 
son, in the case of noisy data, one has to look for a regularized approxima- 
tion of the C-generahzed solution, which can be obtained by minimizing the 

functionalt^^l'f^^] 

*a[V1 - \\LV - g\\\ + X\\CV\\\. (9.26) 
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If we denote by Vx the minimum of the functional (9.26) and if we put 

Vx = il^x(s)t + (l>x(s)n (9.27) 

then it is easy to show that, on each regular arc, V'A and (px must be sohitions 
of the system of differential equations 

-h(s) + 2e(s)4>x(s) -f \0(s)f'i^x(s) + 0{s)<t>x(s) - (9.28a) 

-A[(;iA(5) + 20{s)il^x{s) + \e{s)\'cf>x{s) + 0{sU^x(s)] + <l>x{s) = g{s) (9.286) 
plus boundary ^onditions similar to those discussed in the previous case 
(continuity of V'a at the discontinuity points, etc). The determination of 
the parameter A can be performed using one of the methods discussed in 
Section 5. 

In practice, the most economical method for the computation of Vx 
is perhaps the conjugate gradient method. Regularizing properties of this 
methodf^^ 'f^'^ can also be used in order to avoid the minimization of (9.26). 

In the previous treatment we have neglected the errors in the determi- 
nation of the contour which imply an approximate knowledge of the operator 
L (Equation (9.8)). However, if the equation LV = g + dg, where dg \s the 
error on the data, is replaced by the equation (JD + dL)V = g^dg, where dL 
IS the error on the operator, it appears that the two equations are equivalent 
in the sense that only the error on g is different in the two cases (in one case 
it is dg and in the other case it is dg - {dL)V). This point of view assumes 
that the errors in the determination of the contour have been included in the 
errors on the data. 

9.3. Two-dimensional optical flow 

As we already recalled at the beginning of this section, Horn and Schunckf^*^] 
attempted to recover the optical flow in the entire image and not just on 
a one-dimensional contour. Their basic equation is Equation (9.1), which, 
written expHcitly, provides the relationship 

VE-V^-dtE (9.29) 

where VE = {d^E.dyE} is the gradient of the brightness distribution in 
the image, V is the velocity field (optical flow), and dtE is the partial time 
derivative of the brightness. Therefore, a measurement of VE and dtE gives 
the component of V parallel to VE. 

We assume that the brightness distribution E(x,y,t) is defined in a 
bounded region Q whose boundary dU is a contour with an everywhere con- 
tinuous tangent. Furthermore, we will also assume, for simphcity, that VE 
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is never zero in ft and that the level lines of E{x,yJ) have everywhere differ- 
entiable tangent and normal. We denote by fand n the tangent and normal 
to the level line at the point {x,y} 

-IV^l-'jif^) , "-|V^r>(Jf). (0.30) 

Then the velocity field V{x.y) can be everywhere represented as follows 

V(x,y) =v^{x,y)f-{- v^(x,y)ii. (9.31) 

The problem can again be formulated as the inversion of a projection oper- 
ator: taking X = Y = L^(n) © L'^(n) and 

{LV){x,y) = v^{x,y)n. (9.30) 

The data will be given by 

9{^,y) = 9{x,y)n, (9.33) 

where g{x,y) is the measured value of ~dtE/\VE\. Then the set of solutions 
of the equation LV = g is the set of velocity fields 

V(x,y) = 4ix,y)f+g{x,y)n (9.34) 

where is an arbitary function in L'\Q). The generaHzed solution V'^ is 
trivial also in this case, since V^ = g. 

9.4. A C-generalized solution for the two-dimensional optical 
flow 

As in the case of the optical flow along a contour, it is necessary to look 
for C-generalized solutions. The method suggested in [391 Horn and Schunck 
(1981) can be formulated in this framework. 

Introduce the constraint space Z = X ^ X and define an operator 
C : X H-> Z as 

with the associated seminorm 

\\CV\^z = f {d^V ■ d^V + dyV ■ dyV}. (9.36) 

Written in terms of the cartesian components of V this is just the integral 
of the quantity called the measure of the departure from smoothness in the 
velocity flowf^^^ . 

First consider the question of uniqueness. The null space N{C) is the 
set of the constant velocity fields, say V = a, while the null space N{L) is the 
set of the velocity fields which are orthogonal everywhere to n, i.e., V -n = 0. 
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The intersection is the set of constant velocity fields such that a- n = and 
this condition cannot be satisfied by a ^ if the level lines are not parallel 
straight lines everywhere. 

It is easy to verify that conditions (i) - (iii) of Section 3.2. are satis- 
fied and the existence of the solution is guaranteed. It may be interesting 
however to write the Euler equation for the C-generahzed solution. After 
some lengthy but elementary computations, using the orthogonahty rela- 
tions n . d,n = ri • OyU = 0, t- dj= t- dyt = 0, d,n • dj = dyti ■ dy{= we 
obtain 

\\CV\\% = ^ {iV^l' + [An]' + \dynf) \gf]dxdy-^ 
^ {iVvi^ + (15./? + \dyt]') \U'\' + 2 [(.7 . dj-) d^g + {a • dyt) dyg] ^^ + 

+2g [{{• dj~) d,v + {t- dyu) dy4^] \dxdy. (9.37) 

In order to find the Euler equation of the functional (9.37) one has to consider 
a variation of ^', V ^ V' + h and put equal to zero the term of first order 
in h. Then, using the divergence theorem in order to eliminate the partial 
derivatives of h, transforming the fourth term in Equation (9.37) by means 
of the identities t- d^n = -n • dX t' dyU = -n • dyf, and using the fact 
that h is arbitrary, one finds that the unique function which minimizes the 
functional (9.37) is the unique solution of the boundary value problem: 

VV'+(|5./l' + |a,/f)04- 
+2[(a . d^f) d^g + {n ■ dyt) dyg] + (n • At) g = (9.38) 

di!> / dt\ 

where u is the normal to dQ. Notice that this boundary value problem is 
just the extension in the 2-D case of the problem (9.19) with the boundary 
conditions (9.21). The boundary condition (9.39) can be replaced by the 
value of if the tangent velocity can be measured on dfl. 

It is also easy in the present case to verify that if the motion is a pure 
translation (i.e., a constant velocity field), and if the data function is noise- 
free, then the C-generahzed solution coincides with the exact velocity field. 

It is also obvious that in such a case the C-generahzed solution is ill- 
posed and that one must introduce regularized approximations. These can 
be obtained by minimizing the analogue of the functional (9.26), and this is 
precisely the method used in [39]. 
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10. Surface reconstruction 



Most algorithms able to recover depth from pairs of stereo images provide 
sparse depth values; that is, depth is obtained only for special points in 
the viewed scene. Since a global description of the 3D structure of the 
viewed scene is desirable, it is useful to consider the problem of recovering 
a mathematical representation of a visible surface f{x,y) from the sparse 
datat^«''[^i]. 



10.1. Surface interpolation 

The original data are a finite set of depth values -, =. /(^,,/y,), / = 1 , ,, 

(which are assumed to be exact; that is, noise-free) and the problem is the 
recovery of a smooth function f{x,y) interpolating z, at {x,,y,) ^ f, con- 
tained in n. Grimsont^*^ 't4ii proposed to find / such that it minimizes the 
seminorm 






dxdy. (10.1) 



Uniqueness of solution is guaranteed by the existence of at least four non- 
coplanar points Zi = /(.r,, t/,) f'*ol,[4i] 

This procedure can be seen as an application of generalized inverses in 
the case of discrete data (see Section 3.4): in this case, uniqueness of the 
solution is guaranteed when the intersection of the null space of C,{N{C)) 
and the null space of L{y{L)) is empty, where L is the operator defined in 
Section 3.4. 

The null space N{C) is composed of the set of functions f{x,y) = ax + 
by + c with a, 6, c constants. These functions consist of all planar surfaces 
defined in Q. The null space N{L) has been defined in Section 3.4 and 
consists of the set of functions such that f(x,,y,) = for i =: l,...,n. 
Therefore, it is easy to see that when / ^ 4 and the points /, =:: (.r,,{/,) 
are distinct, the intersection of N{C) and N{L) is empty. In other words, 
uniqueness is guaranteed if there are at least four non-coplanar points, as 
required in [40], [41]. 

10.2. Surface approximation with noisy data 

It is also useful to consider the case in which the data are noisy; that is, when 
the original data are g, = f{t,) 4- £,, i = l,...,N and s, is additive'noise. 



40 



In this case, it is reasonable to look for a solution close to the original data 

n hilt alcio «Tnr,r.fli[3!^]440],[4ll,[42],[43lJ44l rrii • , , 

^,, out also smootUL i ^ i-^ il j xj^jg approach can be seen as an 

appUcation of regularization theory. In Part ne we have already shown that 
interpolation is an ill-posed problem which can be solved by the use of a gen- 
eraHzed inverse. We will present now an approach to interpolation directly 
originating from regularization theoryf^'^J'f^^l, which clarifies the relationship 
between splines, regularization theory, and gives a different framework to the 
results on visual interpolationt^'^J'f'*il'[42],[43],[44],[45] 

We can consider the case in which we want to estimate a smooth function 
f{-t)J e Q C R\ given a finite number of observations of linear functional 
of /. In the case of spatial interpolation, our functionals are: 

9i = F,{f) + ei = /(//)-f e, i = 1,...,/?, (10.2) 

where t, is additive noise. A regularized estimate fn,\ is obtained by solving 
the minimization problem 

E(/(^^)-^«j ^^Jmif). (10.3) 

in which J^() is a seminorm in Hm (H^ is a reproducing kernel Hilbert 
space of functions defined in 17) defined by 

(here m indexes the highest square integrable derivative) and A controls the 
tradeoff between the degree of approximation of the solution to the data and 
the smoothness of solution. The value of A can be computed by the method 
of generalized cross-validationf2il't22] j£ ^^^ ^ 2 we have the functional (10.1). 
The solution of this minimization problem is one of the "thin plate spHnes," 
so called because J-iif) is the bending energy of a thin plate. 

In [58] it was shown that a unique solution exists for any A > provided: 

(1) m >1; 

(2) n -A/-('"+^); 

(3) the "design'^ /i, ...,/„ is unisolvent, that is if {(t>^)'l^^ is a basis for the 
M dimensional space of polynomials of total degree m - 1 or less, then 
S"=i "..^^(/t) = 0(2 = 1, . . . ,n) implies that the a^ = 0. 

If m =^ 2, then we need at least three points which do not He on the 
same straight Hne (to satisfy the requirement of a unisolvent design), which 



41 

is the same requirement as found in [40] and [41]. JMoreover, the solution has 
an explicit representation^^^^ as: 

"( n 

/n,m.A(0 = 5^Cj£'„,(/,/^) + ^(/^(;6^(#), (10.5) 

where 

E,n{sJ) = 0J^s -t\'^\og\s - t\ 
with 

'^ = (-ci,fyi) ^ = (.r2,^2) \s -f\ = \/(xi - .C2)'^ + (yi - y-iY 
and 

^„, = l/2^"'-i7r[(m-l)!]'. (io.7) 

The coefficients c = (c, , . . . , c„) and (i = ((fj , . . . , cf„) are determined by 
the solution of the algebraic linear system: 

T^C = (10.8) 

where K is the n x n matrix with Kjk = Em(tj,fk.). p -= n\,T is the n x m 
matrix with T^, = 0^(/,) and ^f = (g^ , . . . ,^„). 

10.3. Surface interpolation on a regular grid 

While surface interpolation from sparse data requires an arbitrary grid of 
knots, other problems of machine vision require the approximation of a 3D 
surface through points given on a rectangular grid. For example, when a 
smooth function / interpolating intensity values on the regular grid of a 
CCD camera is regularized, it is possible to use doubly cubic splines or a 
tensor product of sphnes, giving an interpolating function that minimizes 



j j [d'fldx'diffdxdy. (10.9 



In this case different kinds of doubly cubic spHnes can be used, according to 
the available data^^*^'. The algorithms are then convolution algorithms (see 
section 8.4). 
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11. Shape from shading 

It is a common experience to notice our ability to recover the shape of an 
object from its shading. Convexity or concavity of viewed objects are easily 
understood by looking at the profile of radiating light. Here we have another 
classical problem of early vision, "shape from shading," which has stimulated 
elegant mathematical approaches. The problem of shape from shading was 
initally formulated in [61] and [62] as the solution of five ordinary differential 
equations called the characteristic strip equations. Of considerable use in 
this problem has been the reflectance map R(p,q) [^^I't^^] ^jji^-jj specifies the 
radiance of a surface patch as a function of its orientation, determined by 
the pair {p,q). If z{x,y) is the surface of the object, p and q are defined as 

dz ^ dz 

and the unit normal n to the surface is 

" " /-, , ■> . {-P^-ga}. (11.2) 

The reflectance map can be computed from the bidirectional reflectance- 
distribution function and the light source arrangementt^^l. 

Formally, given an image E(.r,t/) and a reflectance map R{p,q), the 
shape from shading problem may be regarded as the recovery of a smooth 
surface z(x,y) satisfying the image irradiance equation 

over some domain 1) of the image. Since there are two unknown functions 
{p and q) and only one equation the solution is not unique and the prob- 
lem is underconstrained (and ill-posed). Uniqueness of the solution can be 
recovered by the use of photometric stereo, which takes multiple images of 
the same scene from the same position with diflFerent illuminationf^^J . In this 
approach, several equations of the type of Equation (11.3) are available, with 
different reflectance maps since the illumination source is different. Three 
different light sources can be used to obtain a unique solution. 

If only one source of illumination is available, uniqueness can be restored 
by variational techniques similar to those previously seen. Assuming that the 
object has a Lambertian surface and is illuminated by a planar wave of light 
(and the unit vector s = (51,^2,53) points to the Hght source), then the 
Lambertian reflectance map becomes 

R{p,q) =n- s. (11.4) 



p^ + q^ 
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If, instead of using the pair {/>, g}, the new variables {f.g} are introduced 

J = - ., 9 = , (11.5) 

the reflectance map becomes 

4/ 4^ 

ai-i* (11.6) 



4-(/'^+5'')' 4-(/2+^2) 
The problem of shape from shachng can be formulated either using the un- 
known n or the pair {p,q\ or {f,g}. 

11.1. The variational approach to shape from shading 

When the unknown n is used, the variational approach is to find n{x,y) such 
that it minimizes 

/ (E{x,y)- n-sydxdy^\ I {nli-nl)dx dy, (11.7) 

with the constraint | :ri|| =1. In this case, the variational problem is quadratic 
in the unknown n, but the constraint ||r7|| = 1 is unusual. 

When the pair {f^g} is used, we seek functions / and g minimizing: 

/ \{E{x,y)-R{f,g)\'dxdyi-\ f {fl i- Py ^ gl + gl)dx dy, (11.8) 

with R(f,g) given by Equation (11.6). The variational problem is not 
quadratic in the unknown {f,g} and the results of nonlinear inverse problems 
have to be used. 



11.2. Regularization of shape from shading 

We give an apphcation of the result stated in Section 6 by formulating the 
problem in terms of the pair {p,q}. We define the space X as the direct sum 
L'^Q) © I^(n), i.e., u is a pair {p,q} of square integrable functions: 



I 2 



X 



/ P^{x,y)dxdy i- / q^(x,y)dxdy. (11.9^ 

Jn Ju 



Let the space 1' be also a space of square integrable functions (we now call 
g{x,y) the image E{x,y)), and from (11.2), (11.4) we define a nonlinear 
operator .4 : A" — > 1' as follows: 

i A \i \ -^3 ~ p^\ — q^2 

iAu){x,y) = ' ^ \ (11.10) 

V 1 + p^ + ^2 
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Because n and 5* are unit vectors, it is obvious that j(.4«)(.r{y)| :-. 1 for any 
{.r, (/} e Q. It follows that the domain of .4 is X and that the range of .4 is 
contained in the set oi g{x,ij) such that \g(x,y)\ < 1 in Q. Furthermore, it 
is not difficult to prove that the operator .4 is continuous everywhere, i.e., if 
u is any element of A" and if {tt„} is a sequence convergent to u, then Aun 
converges to An. Indeed, using the inequalities 

1^3 -psi -qs2\ < Vl+p2 +q2^^i +p2 4_^2 > i^ (11.11) 

it follows that 

\AU ' AlLnl < \si\\p-p„\ + \s2\\q-qn\ + | \/l + p^ + ^^ _ ^^ + p2 _|. ^2 | _ 

(11.12) 
Then, using the inequality (^i ^ ... + 9^)2 <^ n{qj + ... + (?,;) (with n = 2,3), 
we get 

\AU -AUrrl'^ < {\P~Pn\'^ + l^ - qj^) ■ (11.13) 

By integrating over Q we get the continuity of the operator .4. 
Finally, we consider the constraint operator C defined by 

l'^"iiz = / [Co{p' + q') + C,(pI ^pl^ql^ ql)]dxdy (11.14) 

where Co could take the value Co = and give the stabiHzer used by Ikeuchi 
and Horn. 

We can seek a solution to the problem of shape from shading by mini- 
mizing the functional 



/ 

Jn 



\{Au)(x,ij) - g(x,y)\2dxdy ^ X\\Cu\\l, (11.15) 

'n 

where the first term in Equation (11.15) is Equation (11.10) and the con- 
straint operator C is defined in Equation (11.14). Because the operator .4 
is continuous and the constraint operator has a compact inverse, the results 
presented in Section 6 indicate the existence of at least a local minimum of 
the functional (11.15). Furthermore, ii g E R(A), the ux converges to an 
exact solution when A ^— > 0. 
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12. Discussion 



The review of early vision presented in Part Two shows that certain regular- 
ization techniques can be useful for a correct and sound "solution" of many 
vision problems. The key idea of all regularization techniques is to introduce 
a priori knowledge — or constraints — which have to satisfy the solution. 
Therefore we will have different solutions according to the assumptions we 
have made, that is our a priori knowledge about the world. 

Physical plausibility of the solution is the most important criterion in se- 
lecting a good solution. The decision regarding the choice of the appropriate 
stabilizing functional cannot be made judiciously from purely mathemati- 
cal considerations. A physical analysis of the problem and of its generic 
constraints plays the main role. Standard regularization theory provides 
a framework within which one has to seek constraints that are rooted in 
the physics of the visual world. Standard regularization, however, offers a 
restricted universe of possible constraints since only certain a priori assump- 
tions can be translated into the language of Tikhonov stabiHzers. 

In our example of the computation of motion, the constraint of smooth- 
ness is justified by the observation that the projection of three-dimensional 
objects in motion onto the image plane tends, in a probabiUstic sense, to yield 
smoother velocity fields t^'l. In the case of edge detection the constraint of 
a small norm for the derivative of image intensity is directly justified by the 
bandHmiting properties of the optics. In the case of motion, however, and 
more dramatically in the case of surface reconstruction, the constraint of 
smoothness is not always correct. This suggests - as we will discuss later - 
that more general stabiHzing functionals are needed to deal with the general 
problem of discontinuities. 

A method for checking physical plausibiUty of a variational principle is 
to use the Euler- Lagrange equation associated with the variational problem. 
In the computation of optical flow, the following sufficient and necessary 
condition has been obtained t^^^ (see also Section 9) for the solution of the 
variational principle Equation (9.26), to be the correct physical solution: 

''^ = '^ (12.1) 

where / is the tangent vector to the contour and V is the true velocity 
field. The equation is satisfied by a uniform translation or expansion and by 
rotation only if the contour is polygonal. Therefore, the smoothness principle 
will give correct results when (a) motion can be approximated locally by pure 
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translation, rotation or expansion, or (b) objects have images consisting of 
connected straight hnes. In other situations, the smoothness principle will 
not yield the correct velocity field, but may yield one that is qualitatively 
similar and close to human perception^^^'t^^L In the corresponding case for 
edge detection (intended as numerical cUfFerentiation), the solution is correct 
if and only if the intensity profile is a polynomial spHne of appropriate degree. 

From a more biological point of view, a careful comparison of the vari- 
ous "regidarization " solutions with human perception promises to be a very 
interesting area of research, as suggested by Hildreth's work on the compu- 
tation of motion. For some classes of motions and contours, the solution 
of Equations (9.19) and (9.22) is not the physically correct velocity field. 
In these cases, however, the human visual system also appears to derive a 
similar, incorrect velocity field^^^^'f^^l. 

It may be useful to remember again that Tikhonov stabilizers do not rep- 
resent the only way to regularize ill-posed problems. Different or additional 
constraints such as shape (monotonicity, convexity) have been proposed. 

The most obvious way to solve inverse ill-posed problems without requir- 
ing smoothness of the solution is to use Markov Random Fields as proposed 
in [36] . In this approach, discontinuities can be preserved introducing ap- 
propriate fine processest^^Jf^^J and appropriate potential terms. A possible 
problem in this approach is the critical dependency of the solution on the 
parameters coupUng the different Markov Random Fields. 

12.1. Stereo matching 

Not all inverse problems of early vision can be solved using the regularizing 
techniques introduced in Part One. For example, stereopsis, which is the 
process that computes depth from two images of the same scene obtained by 
two eyes or cameras, appears as an inverse problem that may be approached 
with standard regularization techniques. It turns out that tliis is, however, 
quite difficult. The critical problem in stereopsis is the correspondence prob- 
lem, that is, the matching of corresponding features in the two images. Let 
us consider the 1-D matching problem, by considering the intensity profile 
— or some corresponding feature map — on conjugated epipolar Hnesf^'L 
In this case, the obvious way to match the right image R{x) with the left 
one L{x) is to find the disparity d{x) such that the two intensity profiles 
L{x) and R(x -f f/(.r)) are as close as possible. We can formaHze this in the 
following way: let us define an operator Pr that depends on the image as 

PRf{x)^R{x^f{x)). (12.2) 



47 

The disparity function that we want could be seen as the solution to the 
inverse problem: 

L(x) = PRd(x). (12.3) 

The operator in Equation (12.3) which has to be inverted depends on the 
data and is not known a priori. This class of problems is not covered by 
the available mathematical results. We could still try to determine d{.v) by 
minimizing 

\\L{x) - R[x + d(x))\\. (12.4) 

A sufficient condition for the solution of (12.4) to be unique is that L(x) 
and R{x) are strictly monotonic functions of x. This is clearly a very re- 
strictive condition, almost never satisfied by real images. In general, the 
problem admits many solutions unless constraints are imposed on d(x). If 
we use constraints of the Tikhonov type, we can look for a solution d[x) that 
minimizes 

Ui-^) - R(x ^ d(x))^ +A||/(^)||. (12.5) 

The second term in (12.5) is the disparity gradient, which is thus introduced 
as a direct consequence of regularization methods. 

One important property of the disparity is that d(x) can be discontin- 
uous. Furthermore, there are often occlusions, that is regions in one image 
that do not correspond to any part in the other image. In this case, d{x) is 
not defined. 

Because of the presence of occlusions and discontinuities in the dispar- 
ity, Equation (12.5) does not provide a physically plausible solution. Equa- 
tion (12.5) requires d(x) to be continuous and diff'erentiable. Equation (12.5) 
is, however, valid if the disparity gradient is strictly less than 2 ( Julesz' def- 
inition): in this case there are no occlusions and Equation (12.5) provides a 
physically plausible solution. 

Another problem with Equation (12.5) is that in many instances match- 
ing is not performed between the intensity profiles in the two images, but 
between features maps. In this case, L(x) and R{x) are not continuous 
functions of x. 



12.2. Pseudoinverses versus regularized solutions 

It may be useful to summarize some mathematical conclusions on the re- 
lationship between pseudoinverses and regularized solutions that may be 
relevant in early vision. 
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(1) When we are dealing with operators defined on Hilbert spaces X,Y 
of the type 

L:Xy-^Y, (12.6) 

we have three main cases: 

(a) if L is injective, hnear, and continuous, and the range of L is given 
by R(L) = F, the inverse problem is trivially well-posed because the 
inverse operator is continuous; 

(b) if L is not injective, and R{L) is closed, then by the method of pseu- 
dosolutions, the inverse problem becomes well-posed in the sense that 
the generalized inverse is continuous; 

(c) a R{L) is not closed, the use of pseudosolutions by itself does not guar- 
antee the existence and continuity of the inverse solution in the case of 
noisy data; then, we have to apply regularization methods. This case is 
the normal one in early vision problems, because noise is always present 
in the data. 

(2) When we are considering operators defined on finite dimensional 
spaces i?" and i?'", of the type 

L : R" ^ R'"" (12.7) 

we have again three main cases: 

(a) If p denotes the rank of the matrix associated with the operator L and 
p= a ^ m, then L is injective and R{L) = R"\ A solution always exist 
and is unique and the inverse problem is well defined. 

(b) If p < ju then uniqueness does not hold, but it can be restored by 
considering the generalized Moore- Penrose inverse. 

(c) If p < m, then existence does not hold for arbitrary data but it can be 
restored again by considering the Moore- Penrose inverse. 

For operators in finite dimensional spaces, the inverse or generalized in- 
verse is always continuous. Therefore with finite dimensional spaces, the 
use of generalized Moore-Penrose inverses is sufficient to guarantee well- 
posedness and the use of regularization methods is not strictly required. It 
is useful to remember again that well-posed problems can be ill-conditioned, 
and in such a case it is necessary to use regularization methods as in the case 
of ill-posed problems. This is the case when the input data are very noisy 
and when differentiation on the input data is required, as in the recovery of 
optical flow or edge detection. Finally, it is useful to observe that: 

(3) The distinction between interpolation and approximation of discrete 
data is associated with the use of pseudoinverse solutions or regidarization 
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methods. Regularization methods are intrinsically approximating solutions, 
while pseudosolutions can be seen as interpolating solutions through the 
original data. 

12.3. Learning regularization algorithms from examples 



12.4. Limits of regularization methods in vision 

We believe that regularizing methods provide a useful and mathematically 
sound framework for many problems in early vision. Not all problems in 
early vision, however, can be easily treated in this framework (e.g. stereo 
matching). Moreover it is often desirable to have solutions that preserve 
discontinuities. In this case Markov Random Fields formulations'^^J are hkely 
to represent more powerfid tools. 

Beyond early vision, efficient vision systems are hkely to use more com- 
plex and elaborate a priori information and to be equipped with reasoning 
capabiHties which encompass early vision and regularization methods. It is 
possible, though as yet unclear, that between early vision and high vision 
different sources of low-level informations are integrated into a unified repre- 
sentation of surface properties, such as the 2^D sketch. At this point MRF 
models could be quite useful in providing a flexible (though complex) tool 
for sophisticated integration. 
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13. Appendices 
13.1. Appendix A. 

For the convenience of the reader, we summarize in this appendix some basic 
ideas of functional analysis which are used in the paper. 

All the questions of existence, uniqueness, and continuity of the solution 
have a precise meaning (and a precise answer, when the answer is known) if 
we carefully define the sets X and Y to which the functions u, g (Equation 
(3.1)) belong. In particular, continuity requires that we define what we mean 
by the "vicinity" of two functions in X or Y; i.e., we must introduce a metric 
in both spaces. 

Among various possible choices of metrics, the one corresponding to a 
Hilbert space is the most simple and interesting, since a Hilbert space is the 
space most similar to the usual Euclidean space. 

A Hilbert space X is a finear space of functions, satisfying the following 
conditions: 

(a) For any pair of functions u,v t A^ (with * being complex conjugate, a 
complex (real) valued function {u,v)x^ called the scalar product in X, 
is defined such that: 

(i) («, u)x > 0, = 0, if and only if u = 0; 

(ii) {a,v)x = ii',ii)\; 

(iii) {Xu + f.iL\z)x =^ \{ii,z)x + u{l\z)x for arbitrary complex (real) 
numbers A, yu; 

(b) A" is complete; i.e., the Cauchy criterion for the convergence of sequences 
holds true; 

(c) A^ is separable; i.e., there exists a countable orthonormal basis. 

The classical example of a Hilbert space is provided by any space of 
square integrable functions (I^-spaces), the scalar product being defined by 

{a,v)x - / ii{x)v^x)dx. (Al) 

A norm can be introduced in A" by means of the relation 

Ii II / \i/2 

||«ILy = {u,ii)x ^ (.42) 

and it satisfies the properties: 

(i') \\u\\x > 0,= if and only if ?/ = 0; 



p{u)= ( / \u'{x)rdx) ' . (.441 
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(ii') i|Au|l Y = |A|||f/|I.Y for any complex (real) A; 

(iii') \\u 4- t'||,Y < I!"I!a' ^ i!^'!Lv (triangle inequality). 

Then the distance q[u,v) from u to v is defined by 

0(11, v) = \\u - v\\x. (.43) 

A real- valued functional p{il], defined on A", is Ccdled a seminorm on A", if it 
satisfies properties (ii') and (iii') of the norm. Then it follows that />(0) = 
and p[u) > but the condition p(u) = does not necessarily imply u = 0. 
The set of the functions suck that p{u) == is a Hnear subspace called the 
null space of p, i.e., it contains the zero of A", and, if it contains u,v, then it 
contains also \u-j-/.iv, for any complex (real) A, /i. An example of a seminorm, 
which is not a norm, is the following: 

( / I w'(a!)| dx 1 

The null space oi p(u) contains all the constant functions. 

An operator from a Hilbert space A" into a Hilbert space Y is defined 
by a mapping which transforms functions of a subset of X (the domain of 
the operator, denoted by D{L)) into functions of a subset of Y (the range 
of the operator, denoted by R(L)). When the domain is a Hnear subspace 
and the mapping is linear, the operator is called Hnear. We use the notation 
L : X —^ Y ior denoting a Hnear operator from A" into Y. 

If D(X) = A" and if, for any sequence «„ converging to it, the sequence 
Lun converges to L«, then I is a continuous operator. A Hnear operator L 
is continuous if and only if it is bounded, i.e., there exists a constant C such 
that, for any u E A" 

\\Lu\\y' < C\\a\\x. (.45) 

The quantity 

,,_, sup ||Xit||y' 

a t A \\u\\x 
is called the norm of the operator L. 

A Hnear continuous operator L : X -> Y always admits an adjoint 
operator L* : 1' — > X, defined by 

{Lii,r)Y ={u.,L*v)x (.47) 

for any 11 e X, v t 1'. The operator L* is also Hnear and bounded and 

I I F 4c I I I I T~ I I 

Using the definition of the adjoint operator, we can write Equation ( A6) 
in the following form 

sup {L*Lu,u)x]'^' 



L 



G X ( u , li ) \ 



.48] 
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sic 
therefore, [|Zj| is the square root of the supremura of the spectrum of X L. 

Here we have used a result which is an extension of the well-known variational 

property of the maximum eigenvalue of a matrix. 

A particular class of linear operators is provided by the hnear function- 
als, which correspond to the case where Y is the space of the real (complex) 
numbers. Therefore, a functional is an operator which associates numbers 
to elements of A". Linear continuous functionals are characterized by the 
Representation Theorem of F. Riesz: let F[u) be a continuous functional on 
X, then there exists a unique function (p ^ X such that 

F(u) ^{u,(p)x (.49) 

for any u G A". 

A Hnear operator L : X ^ Y \s said to be compact (or also completely 
continuous) when it is bounded and transforms any bounded set of A" into a 
precompact set of 1' (a precompact set is a set whose closure is compact; i.e., 
it has the Bolzano- Weierstrass property). Compact operators are interesting 
since they are the most similar to matrices. 

In order to establish the compactness of an operator, one needs com- 
pactness criteria in functional spaces. The basic result is the Ascoli-Arzela 
Theorem, whose proof is the paradigm of all the proofs of compactness. 
Ascoh-Arzela's theorem states that a sequence of continuous functions Un{x) 
is precompact if: 

(i) The functions «n(.r) are uniformly bounded; i.e., there exists a constant 
C such that 

WnM\ < C (.410) 

for any n and any x. 

(ii) The functions u„(.c) are uniformly continuous; i.e., for any t :• there 
exists 6 > such that 

\Un{x) - U„{x')\ < £ {All) 

for any x, x' with \x - x' \ S S and any n. 

A very simple example of a compact operator in a I^ -space is provided 
by an integral operator 

{Lii){x) = J K{x,y)u{y)dy, (.412) 

whose kernel K(x,y) is square integrable 

dx dy\K{x,y)\'^ < +oo. (^13) 
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The most striking property of compact operators in Hilbert spaces is 
that they have a singular value decomposition (SVD) similar to that of a 
matrix. The singular system of a compact operator is defined as the set of 
the solutions of the coupled equations 

Luf, = akVf, , L*Vk- = akUk, (-^14) 

where the aj, (the singular values) are positive numbers and the Uk,Vf, (the 
singular functions) are functions in X and Y respectively. 

When L is compact, its singular system always exists and has the fol- 
lowing properties: the a^ have finite multipHcity and tend to zero when 
A: -> oo (we exclude here the case of a finite rank operator); the u^ form an 
orthonormal basis in the orthogonal complement of iV(Z) and the v/, form an 
orthonormal basis in the orthogonal complement of N{L*), i.e. the closure 
o{ R{L). 

It is easy to verify that a function g G Y is in the range of L if and only 
if (Picard conditions) 



g^NiLn- , 5] l(^i!:|lii! < +OC 



A. ^^ 



(-415) 



Since a^ ^ when L is not a finite rank operator, it is clear that an arbitrary 
function orthogonal to N{L*) does not always satisfy the second condition 
in Equation (A15) and therefore R(L) is not closed. 

We conclude that the problem (3.1) with L compact is ih-posed. 

On the other hand, the range of a finite rank compact operator is closed 
since its dimension is finite. 

Another example of continuous operators with closed range is provided 
by the projection operators; i.e., Hnear operators P such that: 

(i) P* = P; 
(ii) P2 = p 

It follows that liPlI =. 1. Furthermore, N(P) is the set of all the functions 
u = (I- P)u, where v is an arbitrary function of X and R(P) is the set of all 
the functions u such that u = Pu. Therefore R{P) is closed. A projection 
operator is compact if and only if the dimension of R{P) is finite. 
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13.2. Appendix B 

111 this Appendix we will outline the proof of the main result stated in Section 
3.3; namely, that if the constraint operator C satisfies conditions (i) - (iii), 
there exists a unique C-generalized solution aj for any g such that Pg t 
LD{C). However, we first show that conditions (ii) and (iii) are satisfied by 
seminorms defined in terms of differential operators. 

Consider in i'^(0,l) the following seniinorm: 

P(u) = (^j \u^'\x)\'dx^ , (51) 

where t/^^^ = d^u/dx^. The domain of the operator Cu = u^^^ is the set 
of all functions u which have square integrable derivatives at least up to 
order k. Therefore, C is not everywhere defined in 1^(0,1) and is not con- 
tinuous (bounded). However, a diff"erential operator such as C is a typical 
example of a closed operator, i.e., of an operator satisfying the following 
conditions^^^]; if {{*„} C D{C) is a sequence convergent to u and such that 
{Cu„} is a sequence convergent to r G Z, then u e D{C) and Cu = v. 
Furthermore, in our specific case, R{C) = 1^(0,1) since, given an arbitrary 
function v G L^(0,1), there always exists a function u G D{C) such that 
u^ '' = 1'. Therefore, C satisfies condition (ii). This condition impHes that 
C has a bounded generalized inverse C^ since C+ is just the inverse of the 
restriction of C to N{C)^ (the inverse of a closed operator is also closed and 
an everywhere- defined closed operator is bounded). 

As concerns Condition (ii), notice that N(C) is the set of all the polyno- 
mials of degree <k-l. Therefore, N(C) is a ^^ dimensional closed subspace. 
It follows that, whenever I is a Hnear, continuous operator defined on A" 
with its range in an arbitrary Hilbert space F, LN(C) is a ^--dimensional 
closed subspace in Y . 

Finally, as concerns Condition (i), it is satisfied whenever N{L) contains 
no polynomials of degree v k ~1. 

The proof of the result stated in Section 3.3 is based on the fact that 
conditions (i) - (iii) imply the following one: there exists a constant fS'^ such 
that, for any {/ G D{C) 

\\Lui\l + \\Cu\\l >f3'^\\a\\']^. (B2) 

This condition is called by Morozov the completion condition^'^^l Suppose 
we define on D{C) the scalar product 

{u,v)o = {Lu,Lv)y + (Cu,Cv)z. {BZ) 
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Then from Condition (i) it follows that ||t/.||o = implies ||«|| = 0, and 
D{C) is complete in the topology induced by this scalar product, i.e., D(C) 
becomes a Hilbert space. Condition (B2) indeed impHes that any Caucliy 
sequence in the topology of D{C) is also a Cauchy sequence in X. 

In order to prove inequaHty (B2), we introduce the space \V = Y (B Z 
and define the operator B : X ^ W as follows: 

Bu = {Aa.Cii}, ueD{C). (B4) 

Then, since B is closed, inequahty (B2) is equivalent to stating that B has 
an inverse B"^ and that R{B) is closed in W. 

The existence of j5~^ is an easy consequence of Condition (i). In order 
to prove that R{B) is closed, let {un} C D(C) be a sequence such that 
{Bun} is convergent in W; we must prove that the hmit belongs to R(B). 
Since {5ti„} is convergent in W, it follows that {.4a„} is convergent in Y 
and {Citrr} is convergent in Z. Put w„ = u\^^ + ui!\ with u^^^ G N{C) and 
Un G iV(C)^. As already remarked, the restriction of C to N(C)^ has a 
bounded inverse and therefore there exists a constant 7^ such that 

l|C«JI^ = !|Cu«„''|li>7^II«i."|l|. (B5) 

This impHes that {a„ } is a Cauchy sequence and therefore it is convergent. 
Let if(^^ be the hmit. Since the operator C is closed, u^^^ 6 D{C) and 
{Cun} converges to Cu^^\ Now we have Lun = Lu\^^ + Lui^^ and both 
{Lun} and {Lun } are convergent. It follows that {Lun^} is also convergent, 
and, thanks to the closure of LN{C), there exists u^^^ 6 N{C) such that 
Lun converges to Lu^^\ By combining all the results we have that {5ti„} 
converges to B{u^^^ + u^^^) and therefore R(B) is closed. 

Now, starting from the completion condition (B2), the proof of the ex- 
istence of the C-generaHzed solution for any g such that Pg G LD{C) can be 
done as int^'l . The proof is just an easy extension of the proof of the classical 
result that any closed and convex set has a unique element of minimal norm. 
Notice that in [27], the C-generahzed solution is called the solution of the 
basic problem. 

When the operator C has a bounded inverse C^^ conditions (i) - (iii) 
are obviously satisfied. In such a case, the C-generalized L+ is given in [3]: 

L:^ = C-\LC-')\ (B6) 

where (IC ^^ is the generafized inverse of the operator LC^^ : Z -^ Y. 

An example of an operator C satisfying this assumption is the following: 
takeX == 1^(0,1) and Z = 1^(0, 1) ^^(o, i) and define C by 

Cu = {u,u} i^Bl) 
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with domain the set of the absolutely coiitiimous functions with square inte- 
grable first derivative. Then, 



\Cufz= I \u{.v)?dx + / \ii'{x)\'djc, 

Jo Ji) 



BS 
and we have a functional of the type (5.4). 
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13.3. Appendix C 

In this appendix we show that the problem of Hnear interpolation is equiva- 
lent to the computation of a generaHzed solution in a suitable space. 

Let .Y be a space of differentiable functions, defined on the interval [0, 1] 
and having a square integrable first derivative. A' is an Hilbert space if we 
define a scalar product by means of the formula 

{u.v)x = ii(0)i;(0) + / u'{x)v'{x)dx. (CI) 

Jo 

Let X t 10, 1] be a fixed, arbitrary point; then, from the elementary relation 

u{x) = u{0)+ I iL'(x')dx' (C2) 

Jo 
it follows that 

a{x) =(u,Q^)x, (C3) 

where 

Qxi'P') = 1 + min{.c,a-'}. (C4) 

Clearly Q-, e X for any .r, and therefore all the evaluation functional 
(i.e., the functionals which associate to a function u its value in a given point) 
are continuous. 

A Hilbert space of continuous functions having the previous property is 
called a reproducing kernel Hilbert space. The reproducing kernel Q(.r,.r') 
is defined by 

Q{x,x') = Q,(x') = Q,,{x), (C5) 

and its name is due to the relation 

{Q..Q.')x = Q{x,x'l {C6) 

Assume now that a function u e X is specified at the points Xi, X2, . . . ,xp,r 
{xn t [0,1]) and let Qi^g-i, ■ - ■ ,9n be its values. The interpolation problem 
(i.e., find u e X such that u{xn) = ^„ forn = 1, . . . , iV) can be formulated, 
thanks to (C3), as the problem of determining u e X such that 

{i^'QxJ = 9n ; n = l,...,iV (C7) 

and therefore it takes the form (3.12), (3.13). If we recall that the generahzed 
solution is orthogonal to N{L) (Section 3) and that N{L) is the orthogonal 
complement of the subspace spanned by the functions 

(f>n{x) - Q{Xn,x) (C8) 

{L is defined as in Equation (3.14)), we conclude that the generahzed solution 
must be a Hnear combination of the functions </>„ 

iV 

«^(-P) = ^CnQ{Xrr,x). (C9) 

n=l 
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From Equation (C4) it follows that u^(x) is just the Hiiear interpolation of 
the data gn- 

Interpolation by means of sphnes of degree rn ^ 2k - 1 [k > I) can 
be obtained along similar Hues by a suitable definition of the reproducing 
kernel Hilbert space X ^^'^' . Interpolation by means of natural spHnes of 
the same degree^^^l can be formulated as the problem of determining, in the 
same space, a C-generalized solution which minimizes the seminorm (Bl). 
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